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ABSTRACT. Finite-difference and finite-volume formulations are analyzed in order to 
clear up the confusion concerning their application to the numerical solution of conserva- 
tion laws. A new coordinate-free formulation of systems of conservation laws is developed, 
which clearly distinguishes the role of physical vectors from that of algebraic vectors which 
characterize the system. The analysis considers general types of equations — potential, 
Euler, and Navier-Stokes. Three-dimensional unsteady flows with time- varying grids are 
described using a single, consistent nomenclature for both formulations. Grid motion due 
to a non-inertial reference frame as well as flow adaptation is covered. In comparing the 
two formulations, it is found useful to distinguish between differences in numerical meth- 
ods and differences in grid definition. The former plays a role for non-Cartesian grids, and 
results in only cosmetic differences in the manner in which geometric terms are handled. 
The differences in grid definition for the two formulations is found to be more impor- 
tant, since it affects the manner in which boundary conditions, zonal procedures, and 
grid singularities are handled at computational boundaries. The proper interpretation 
of strong and weak conservation-law forms for quasi-one-dimensional and axisymmetric 
flows is brought out. 


INTRODUCTION 

Many current algorithms in computational fluid dynamics are based on the nu- 
merical solution of conservation laws. This choice is motivated by several consid- 
erations, the chief one being the ability to treat flow discontinuities automatically. 
In an earlier paper (Ref. 1), the author treated the differential formulation of the 
conservation equations. This formulation forms the basis for flnite-difference algo- 
rithms, which are historically very old. The more fundamental integral formulation 
is the basis for flnite-volume algorithms. These appeared more recently (probably 
the earliest is found in Ref. 2), and are not as well known. There appears to be some 
confusion and ignorance concerning the relation of these two approaches. The pri- 
mary purpose of this paper is to shed some light on this subject. For this purpose, 
it is useful to distinguish between the method of approach (integral vs. differential) 
and grid deflnition (grid points represent vertices or centers of cells). It turns out 
that the differences in method only play a role for non-Cartesian grids, and manifest 
themselves in the manner in which the geometry of the grid is handled. The effect 
on the numerical results is generally small in most cases. The differences in grid 
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definition are important at computational boundaries, where they determine the 
manner in which boundary conditions, zonal procedures, and grid singularities are 
handled. 

To provide the greatest possible utility, a unified presentation covering all classes 
of equations is given. Both approaches are described using well-defined, consistent 
nomenclature. The general case of three-dimensional unsteady flow with time- 
varying grids is considered, and specialized to other cases when necessary. Only 
those aspects of algorithms that relate to the comparison of the two approaches are 
presented. 

In order to make the discussion more physical and compact, vector notation is 
used throughout. Here the word vector refers to physical vectors (e.g. velocity or 
momentum) as distinguished from algebraic vectors (e.g. the set of conservative 
variables). It is customary tc replace the vector momentum equation by its scalar 
components, so that only algebraic vectors and matrices are encountered. To allow 
for arbitrary scalar decompositions, and to preserve the compactness, the momen- 
tum equation is here treated as a single vector equation. This necessitates the 
introduction of a novel treatment of flux vectors and flux Jacobian matrices solely 
in terms of physical vectors, and independent of a coordinate system. 

The integral formulation of general conservation laws is presented first, and is 
then used to derive the differential formulation. The new vector formulation of 
flux Jacobian matrices follows, with applications to Roe averaging and flux vector 
splitting. The next section presents the discretization of the equations, beginning 
with a careful discussion of grid definition. The general characteristics of stan- 
dard finite-difference methods are described, and some of the inaccuracies arising 
from the treatment ot the grid geometry are enumerated. The discussion of finite- 
volume methods first details various geometric constructions, including a unified 
presentation of different methods to calculate the cell volume. The three main 
types of equations — Euler, Navier-Stokes, and potential are each discussed in turn, 
and specific comparisons of the two approaches are made. Both centered and up- 
wind treatments of inviscid flux terms are considered. The section on moving grids 
presents in a unified way the combined effects of a non-inertial reference frame 
and grid motion with respect to that frame due to flow adaptation. The role of 
the two types of grids in formulating wall boundary conditions, zonal boundaries 
and grid singularities is covered in the next section. The following section discusses 
strong and weak conservation-law forms, and their relation to quasi-one-dimensional 
and axisymmetric flow. The concluding section contains brief discussions of hybrid 
methods and the relation of finite-element and finite-volume methods. 

This work was partially supported by NASA Ames Research Center under grant 
NCC 2-16. 
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Formulation of Conservation Laws 


Integral Formulation. 


Let the position vector r of a point in space and time t be defined with respect 
to an inertial reference frame. Since physical conservation only has meaning over 
a finite region of space and a finite interval of time, we divide the fiow region into 
contiguous cells which can vary with time. The general form of a conservation law 
is 


j UdV- j 

V(ta) V(ti) 


UdV + 



n-FdS dt 



PdVdt, 


( 1 ) 


where V[t) is the cell volume, ndS{t) is a vector element of surface area with 
outward normal n, C/ is a conservative variable per unit volume, F is the fiux of U 
per unit area per unit time, and P is the rate of production of U per unit volume per 
unit time. If U and P are scalars, then F is a vector; while if U and P are vectors, 
F is a tensor. An example of P is the rate of production of a chemical species. Let 
u and v(t) be the fiuid velocity and surface element velocity, respectively. The fiux 
F can then be written as 

F = (u-v)U + G, (2) 


where the first term is the convection of U relative to the surface element, and G 
stands for the non-convective part of the fiux. It is often convenient to define the 
cell geometry with respect to a non- inertial reference frame. In this case, let ro(t), 
Vo(t),and G(t) be the position vector of the origin, velocity, and angular velocity of 
the non-inertial frame relative to the inertial frame. Then v can be written as 


V=Vr+Vc, (3) 

where 

Vr = vo(t) + n(f) X [r - ro(0], (4) 

and Vc(t) is the surface element velocity relative to the non-inertial frame, but 
expressed in the inertial frame. The velocity Vc can be determined by the motion 
of a well defined surface, or some fiow gradients (adaptive grid). The case u = v 
corresponds to a Lagrangian cell. 

The state of the fiuid is specified in terms of a primaxy or primitive variable Q. 
An example of a fiow governed by a single conservation law is potential flow, for 
which Q is the velocity potential <j> and U is the density p. In general the flow is 
governed by a system of conservation laws, in which case f7, F, P, and Q represent 
a set of variables. The primitive variables Q for the conservation of mass, momen- 
tum, and energy are p, u, and the specific internal energy e. The corresponding 
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conservative variables U are p, the momentum per unit volume m = pu, and the 
total internal energy per unit volume e = p(c+ • u). For nonequilibrium flow, the 
partial density of various species and the energy of internal states of species can be 
additional variables. In some turbulence models, additional variables for turbulent 
quantities may be required. In order to close the system of equations one needs an 
equation of state, and possibly additional relations to define transport coefficients 
and the production terms P. 


For the Euler equations, both F and U are functions of Q, and consequently F 
can be expressed as a function of U. For the Navier-Stokes equations, F depends 
additionally on VQ, where the notation VQ denotes the gradient of Q, its transpose, 
or the divergence of Q. (If thermal radiation is present, F can involve integrals of 
Q, and additional integro-differential equations would be required. This case is not 
discussed in this paper.) For potential flow, both F and U also depend on VQ as 
well as dQfdt in the unstea<ly case. For all flows, an additional explicit dependence 
on r and t comes from the surface element velocity v. We note that the integral 
formulation may require spatial and temporal derivatives of Q to be calculated. 


If we assume all variables are continuous in time, then Eq. (1) reduces to 




/ 


n • F d5 = 


jp,y. 

V 


( 5 ) 


This is the usual statement of a conservation law. For steady-flow calculations, 
the first term is absent. Even when it is employed in time asymptotic marching 
techniques, it need not be treated accurately. This poses less stringent conditions 
on the numerical algorithm. In many applications, Eq. (5) is only satisfied in a 
global sense, over the complete flow region. In this sense it is used as an overall 
check on a numerical method. There are additional global conservation laws that 
can be derived under special assumptions that are sometimes used as constraints 
in an algorithm. Examples are conservation of entropy (inviscid, continuous flows), 
kinetic energy (incompressible, inviscid flow), and vorticity (plane incompressible 
flow). In this paper only local conservation laws that serve to define the basic 
numerical method will be considered. 


Certain geometric identities can be derived as special cases of Eq. (1). 
condition that the cell is closed is given by 


/ 


ndS = 0. 


The 

( 6 ) 


The relative rigid body motion of two frames of reference is expressed by 


/ 


n • Vr 


= 0 . 


( 7 ) 
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Using Eqs. (4) and (6), this can be replaced by 


/ 


r X n dS = 


0 . 


The conservation of volume for a time-varying cell is given by 


V[t2) — V[ti) = j j> n-VcdSdt. 

s{t) 


( 8 ) 


( 9 ) 


An additional geometric constraint is that the sum of the cell volumes must equal 
the total volume of the flow region. 

A main motivation for the conservative formulation is to capture flow discon- 
tinuities in inviscid flow. The jump conditions across such discontinuities are just 
limiting forms of the integral relations. Applying Eq. (5) to a pill box with faces ASi 
and A 52 straddling a surface of discontinuity whose normal in a positive direction 
is n, we obtain in the limit the jump conditions 

j n FdS= j xi‘FdS. (10) 

AjSi A ^2 


Actually, since the geometry of space is continuous, Eq. (10) can be replaced by the 
weaker relation 

n-(Fi-F 2 ) = 0. (11) 

It is possible for algorithms to satisfy the weak form of the jump conditions and 
violate strict conservation. 

Differential Formulation. 

The differential formulation is obtained by applying Eq. (5) to a cell in physi- 
cal space which is the image of the computational cell drj df resulting from the 
coordinate transformation 

r = r(^»/,f,r) 

t = T. 

Let S^drydf be the surface area vector in the positive | direction (with analo- 
gous deflnitions for the other directions), and V d^ dr} df be the volume of the cell. 
Equation (5) then tak.es the form 

{UV)r + (S^ • F)^ -t- (S" • F)„ + (Sf • F)f = PV, (13) 
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where subscripts indicate partial differentiation, 

= r„ X Tj 

S" = Tj X (14) 

X r„, 

and 

F = r^-r„xrj. (15) 

The quantity V is also equal to the Jacobian J or J~^ of the coordinate transfor- 
mation. 

In expressing 17, F, and P as functions of Q and its derivatives, the following 
relations are used: 

V = Tr, (16) 

+ VrjQtj + V^Q^, (17) 

where 

Ve = S«/K, Vtj = S'’/V, Vc = S^/V, (18) 

and 

dQ/dt = Qr-rr‘ VQ. (19) 

From Eq. (19) one can obtain the relations dijd t = —r^ • V^, etc., for the time 
derivatives of the inverse transformation. Equation (13) can also be written in the 
equivalent “Cartesian” form 

Ur + (F()( + (f ’), + (F<), = P, (20) 

where 

U = UV, • F, F" = S" ■ F, Ff = • F, P = PV. 

The geometric identities (6) and (9) become 

(S«)« + (S-), + (S'); = 0, (21) 

and 

Vr = (S^ . r,)e -h (S" • r,)„ + (Sf • r,) j. (22) 

Equation (22) is called the differential statement of the geometric conservation law 
in Ref. 3. Using Eqs. (21) and (22), Eq. (13) can also be written in the quasi- 
conservative form 

17, -h Ve • (F'^ - TrU^) + Vv • (F; - r,t7,) -h Vf • (Fj - r,£7j) = P, (23) 

where 

F' = ut7 + G. 

The above form is called the chain rule conservation law form in Ref. 4. 
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New Coordinate- Free Treatment of Flux Vectors 


We have seen that Eq. (1) can represent a scalar or vector conservation law. Here 
the word vector refers to a physical vector such as position, velocity, or momentum, 
and is represented by a bold faced symbol or an arrow.The related higher order 
quantity is a tensor. It is also convenient to express a set of variables as an algebraic 
vector, and represent it as a column or a row. The related higher order quantity is a 
matrix. Since all calculations must ultimately be done with numbers, and to avoid 
the confusion between the two uses of the term vector, the momentum conservation 
law is normally treated eis three scalar laws for the components of the momentum. A 
compactness and greater physical insight can be obtained by retaining the physical 
vectors as components of the algebraic vectors. The procedure will be illustrated 
for the inviscid flow of a perfect geis. 


Vector Formulation of Flux Jacobian Matrices. 


The set of conservative variables U is given by the column vector 


U = 


P 

m 

e 


(24) 


Let n be the normal in a positive direction to a cell surface or a coordinate surface. 
One can then deflne the normal flux component Fn = n • F, the normal velocity 
components Un = n • u and Vn = n • v, and the normal relative velocity component 
tt' = n • (u — v) = Un — Vn- For inviscid flow, the set of variables Fn is given by the 


column vector 



'm' 


pu' 

Fn = 

P 


mu' + pn 


E 


eu' + pUn 


( 25 ) 


where M, P, and E are the normal flux of mass, momentum, and energy, and p is 
the pressure. 

The vector Fn is a non-linear function of the vector U. In many algorithms one 
therefore employs a linearization in time or space. This requires the calculation of 
dFn in terms of dU. If Vn and n are held flxed, the flrst component of dFn can be 
written as 

d(pu') = —Vn dp + n- dm. (26) 


This can be rewritten in the form of a matrix multiplication as 


d(pu') 


- Vn 


n 0 


dp 

dm , 
de 


(27) 
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where the dot product is implied in multiplying the second element of the row vector 
by the second element of the column vector. Applying the same procedure to the 
other components of dF„, we can define a flux Jacobian matrix operator A satisfying 
dFn = A dU, using the convention that in forming the product of a matrix element 
with a vector element, a dot product is implied if each element is either a physical 
vector or a tensor. For a perfect gas satisfying the equation of state 


p = Kpe, K = 'I — 1, 


the matrix A can then be written as 


A = 


/c6in - «„u un - /cnu + u'l 
(/c6i — H)un Hn — /cu„u 


0 

Kn 

u' + KUn 


( 28 ) 


( 29 ) 


where ^ = 7c + |u • u is the total enthalpy per unit mass, 61 = • u, and I is the 

identity tensor. 

For some algorithms, one requires the diagonalization of A in the form R~^AR = 
A, where A is the eigenvalue matrix, and R and R"^ are the right and left eigenvector 
matrices, respectively. For multiple space dimensions one requires a set of linearly 
independent eigenvectors corresponding to a repeated eigenvalue. Let a, be an 
arbitrary set of spatial basis vectors, and be the set of reciprocal basis vectors 
satisfying a, • a? — where 6,- is the Kronecker delta. One can then define 
a„, = n • a*, b, = n x a^, = n • a-^ , and b-' = n x a-^ . If /? is an arbitrary scalar, 

the three matrices can be written as 


and 



XiSf 

0 

0 


‘u'Sf 

0 

0 

A = 

0 

A2 

0 


0 

u' + c 

0 


0 

0 

A3 


0 

0 

u' — c 


R = 


dfii 

Cn,U + /?bi 
®nt^l / 5 b( • U 


1 1 
u + cn u — cn , 
H + CUn H — CUn 


R-^ 


a{{l - 63) -/ 3 ~^(bJ -u) 

- Un/c) 
libs + Un/c) 


a^^b2U + (3-^h^ -aib2 

-|(62U-n/c) I&2 

-|(62U + n/c) I &2 


( 30 ) 


( 31 ) 


( 32 ) 


where c = is the speed of sound, 62 = and 63 = b^bi. For three 

dimensions, A is a five by five matrix, while R hcis five columns and R~^ has five 
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rows. Since the basis vectors aj are linearly independent, it follws that the three 
eigenvectors corresponding to the repeated eigenvalue Ai are linearly independent. 
One can define functions of the matrix A through 


/(A) = + /(A,)^, + /(As)/”,, (33) 


where the projection operators Pi, P2, and P3 are matrices satisfying 


R-^PiR = 

6/ 0 0 
0 0 0 

, P-‘P2P = 

0 0 o' 
0 1 0 

, R-^PsR = 

0 0 o' 
0 0 0 


0 0 0 


0 0 0 

. 


0 0 1 


(34) 


Examples of /(A) are A, |A|, sgnA = |A|/A, and A^ = (A ± |A|)/2. The formula for 
Pi is 

r 1 — 63 62U ~ ^2 1 


Pi 


u„n - 63U 

u2 - (1 + 63)61 


62UU — nn + I 
(1 + 63 )u - u„n 


— 62U . 

~ ^3 


(35) 


Using the fact that 


Pi + P2 + P3 — / — 


1 

0 

0 


0 

I 

0 


0 

0 , 
1 


(36) 


where I is the unit matrix, one can express the other two projection matrices as 


P, = ^l±A/c+/(l^„'/c)-P,l. (37) 

All the results obtained so far can be generalized to an arbitrary gas law. 

In closing, one should note that the matrices A, P“*, Pi, P2 and P3 are operators 
whose definitions imply that they are operating on some quantity to their right. 
They must therefore be used with some care. As an example, the product R~^R 
gives the identity matrix. On the other hand, the product RR~^ is undefined, and 
does not give the unit matrix I, even when interpreted as a tensor product. 

Roe Averaging. 

Given two states U i and U r associated with the common surface vectors n and 
V, the Roe averaged state U (Ref. 5) is defined by 


FniUR) - F^Ul) = A{U){Ur - Ul). (38) 


This will be an exact relation if the two states correspond to a discontinuity. A 
unique state U can be obtained by invoking the plausible restriction that the average 
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velocity u lies in the plane of Ul and u^. Upon substituting Eqs. (24), (25), and 
(29), and recalling that Ux,, Ur, and n are arbitrary, independent vectors, one 
obtains from the second component of Eq. (38) 

u = aUL + (1 - a)ufi, (39) 

where 

« = v^/(v^+ V^)- (^0) 

Similarly, the third component of Eq. (38) then yields 

H = uHl + (1 - a)HR. (41) 

The sound speed, derived from the total enthalpy via = k{H — |u • u), can be 
written as 

= ^/ca(l - q:)(uh - Ui,) • {ur - Ui,) + ac\ + (1 - oi)cr. (42) 

The above derivation is much more direct than that found in Ref. 5. One notes 
from Eq. (42) that is greater than the weighted average of c\ and c^. It follows 
from Eq. (30) that for either A 2 or A 3 , it is possible that the Roe averaged eigenvalue 
could lie outside the range determined by the two states L and R. In particular, if 
the normal relative speed is very close to sonic for both states, the corresponding 
eigenvalues could both be of one sign, while the Roe averaged eigenvalue could have 
the opposite sign. Some numerical examples illustrating these phenomena are found 
in the appendix. The implications for algorithms based on Roe’s scheme should be 
studied further. 

Flux Vector Splitting. 

In flux-splitting methods the flux vector Fn is split into several parts, each of 
which is approximated numerically according to the sign of a relevant eigenvalue. 
The earliest methods utilize the homogeneity property 


Fn=AU 


(43) 


valid for a gas satisfying 


p = pf{e). (44) 

(It is thus valid for a gas that is thermally perfect, but calorically imperfect.) Using 
Eqs. (33) and (43), one can split Fn as 


Ffi — Fnl + -fn2 + ^n3> 


(45) 
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where 

For a perfect gas one obtains 
Fnl = 


Fni = XiPiU. 


Alp/c 

1 

u 

A a p 

. F^^ = ® 

1 

u ± cn 

7 

Juu 

"3 2 ^ 

H ± CUn 


(46) 


(47) 


This form of flux splitting was first given in Ref. 6 . The earlier Steger and Warming 
splitting (Ref. 7) is of the form 

= F+ + F-, (48) 

where the forward and backward flux vectors and F~ are obtained from Eq. (45) 
by substituting and A^, respectively. If the normal relative speed is sonic or 
supersonic, one obtains 

F+ = F- =0 for u' > c 

F~ = Fn, F^ = Q for u' < -c. 


(49) 


The above split fluxes are not continuously differentiable at the zeros of the eigen- 
values. Van Leer (Ref. 8 ) proposed a new form of splitting which makes each split 
flux continuously differentiable and also degenerate for |u'| < c. The latter con- 
dition produces one zero eigenvalue for the split flux Jacobians, leading to steady 
transonic shock structures with only two interior zones. Equation (49) is still as- 
sumed valid for |u'| > c. For |u'| < c, the requirement of continuous differentiability 
leads to the presence of a factor (u' ± c)^ in the formula for respectively. The 
simplest splitting of the mass flux gives 

M^=±^[u'±cY. (50) 

4c 

The splitting of the momentum flux takes the form 

=M=^[u- i(u'T2c)n]. (51) 

7 

Note that Eqs. (50) and (51) are valid for a general gas law, with 7 = pc^/p. In 
order to satisfy the degeneracy condition, the split energy flux must satisfy 

JE;± = /(M±,P^,n,v). (52) 


One can show that this is only possible for a perfect gas. The resulting expression 
is 


E^=M- 




(53) 
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Discretization of Conservation Laws 


Spatial Discretization. 

A primary grid is defined by choosing a set of discrete grid points and connections 
among them. Normally points are chosen to lie on the flow region boundaries. The 
grid also divides the region into a set of contiguous primary cells, where the grid 
points are now cell vertices, and the connections form the edges of the faces com- 
mon to neighboring cells. If the shapes of the connections are not precisely given, 
then the edges are not defined uniquely. Normally the connections are specified as 
straight lines. Even when the shapes of the connections are given, the shapes of the 
faces are not uniquely defined in three dimensions. For ordered grids the cells are 
quadrilaterals in two dimensions and hexahedrons in three dimensions. Grid points 
may be specified numerically or analytically. In time-dependent algorithms the grid 
points can advance in time. 

A secondary grid can be obtained by determining the centers or centroids of the 
primary cells (in a non-unique way) and connecting them across cell faces. The 
secondary grid points also act as vertices of secondary cells. We thus have two 
interlocking grids, with the cell vertices of one being the cell centers of the other. 
One should note that if the primary grid is not sufficiently smooth, a secondary 
grid with straight line connections may not be possible. One can always define a 
secondary grid with piecewise straight connections by determining the centers of 
the faces of each primary cell, and connecting them to the cell center with straight 
lines. The role of grids in the finite-difference and finite- volume discretization of 
the equations is examined below. The discussion focuses primarily on stationary 
grids at interior points. The case of moving grids as well as considerations at flow 
region boundaries are covered in subsequent sections. 


Finite-Difference Discretization. 


Standard finite-difference methods are based on the discretization of Eq. (20), 
with the geometric quantities included in the definitions of the transformed vari- 
ables. The position vector r and the primary variable Q are defined at the primary 
grid points corresponding to equispaced points in the computational rj, f space. 
The geometric quantities defined by Eqs. (14)-(18) are evaluated at each grid point 
from the r data by central-difference approximations. At interior points, a conser- 
vative difference algorithm is applied in terms of the transformed variables defined 
in Eq. (20). If i, j, k are the integer indices corresponding to the £, t;, f coordinates, 
the spatial difference approximation can always be expressed in the form 




(54) 


where y ^ is any numerical approximation to the flux across the i + |,j,A: 
face of the secondary grid. Similar relations hold for the rj and f derivatives. The 
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numerical flux can be obtained by central differences, or by any of the varieties 
of upwind methods. It can also contain additional artificial terms to stabilize the 
calculation. The time integration of Eq. (20) can be explicit or implicit, may involve 
dimensional splitting, and can involve time linearization to avoid iteration during 
a time step. In the latter instance, the matrix A of Eq. (29) is evaluated at the 
primary grid points. 

The form of Eq. (54) guarantees that when Eq. (20) is summed over all the 
grid points, interior flux terms will cancel. This telescoping property defines con- 

A 

servative differencing. For unsteady flow, the summation implies that Ui,j,k is aji 
approximation to the conserved variable contained in a secondary cell surrounding 
the primary grid point i,j,k. These considerations must be modified for those grid 
points that lie on the flow region boundaries. First of all, the geometric quantities 
must be evaluated using one-sided differences. The flow variables U are calculated 
by applying appropriate boundary conditions. These involve auxiliary equations or 
characteristic relations which determine which quantities are specified and which 
are related to interior quantities in some manner. 

It is seen that the stajidard finite-difference approach combines geometric and flow 
variables and subjects the resultant transformed variables to an algebraic treatment. 
The flow variables are unknown, and can only be determined within some truncation 
errors. On the other hand, for a given coordinate transformation the geometric 
variables obey their own identities which could be satisfied exactly. The combined 
approach can give rise to numerical errors even when the flow is known exactly. 
Some of the sources of inaccuracies are: 

1. The geometric identities (21) and (22) may not be precisely satisfied. Since the 
surface area vectors are linear functions of the position vectors in two dimensions, 
the discrete form of Eq. (21) will have no truncation error in this case when central 
differences are used throughout. However, in three dimensions the area vectors are 
quadratic functions of the position vectors. Since the product of an average is not 
equal to the average of the product, truncation errors will be present in this case. As 
a result, the time integration of a uniform flow will result in numerical oscillations. 
The discretization of Eq. (15) will result in secondary cell volumes that will not sum 
to the total volume of the flow region. This can be a source of error in unsteady flow 
calculations. If the grid is undergoing deformation with time (vc ^ 0), the discrete 
form of Eq. (22) will in general not be satisfied, irrespective of how V and are 
evaluated. This will give rise to errors in a uniform flow even in two dimensions. 
If the numerical flux is obtained by an upwind method, the geometric variables are 
not treated in a fully centered way. This also can produce numerical errors, even 
in two dimensions. 

2. The calculation of derivatives from Eqs. (17)-(19) can give rise to inaccuracies. 
In potential flow, where Q is the velocity potential, the gradient of Q is the fluid 
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velocity and therefore must be determined accurately. If the difference approxima- 
tion for is not consistent with that for based on geometric quantities at the 
primary grid points, the calculation of V(^ from Eq. (17) will result in errors in a 
uniform flow. Similarly, for unsteady potential flow with moving grids, the time 
differencing of Qr and must be consistent in order to avoid further errors from 
the use of Eq. (19). For the Navier-Stokes equations, an inaccuracy of a different 
kind is possible. If a central difference approximation for VQ from Eq. (17) is used 
to calculate the transport terms in the numerical fluxes of Eq. (54), the derivative 
(F^)l at the grid point t,j. A: will depend on the position vectors r at the points 
i + 2, y, k and t — 2,y, k. This lack of compactness puts more stringent requirements 
on the smoothness of the grid in order to achieve a desired accuracy. 

3. The treatment of boundary points can also be a source of inaccuracies. The 
evaluation of geometric quantities by one-sided differences is inconsistent with the 
central difference approximation used at the neighboring interior point, resulting 
in a possible loss of accuracy. If a grid singularity occurs on the boundary, the 
ordinary difference procedure can give a large error. The boundary procedure used 
to evaluate the flow variables at a boundary grid point may not explicitly satisfy 
Eq. (20), which is used at the neighboring interior point. This inconsistency can 
also be a source of error. 

In a practical finite-difference code these numerical errors are handled in two 
ways. Specific corrections can be applied to cancel errors in a uniform flow. The 
standard finite-difference discretization can also be modified to treat the geometric 
terms more precisely. Both of these procedures are discussed in more detail in the 
comparisons with the finite-volume discretization presented below. 

Finite- Volume Geometry. 

As a building block for subsequent calculations, we consider a triangular face 
whose vertices are ri, rj and Vz- A properly oriented surface area vector for that 
face is given by 

Si 23 = <S'i 23 ni 23 = J ndS = -{rz — rj) x (r 3 — ri). (55) 

123 

The formula implies that the edges are straight lines connecting the vertices, and 
that the face is the plane determined by the three vertices. Actually, it follows from 
Eq. (6) that the surface area vector is only a function of the shapes of the edges. 
Thus Eq. (55) is valid for any face with straight line edges. Further reflection reveals 
that there are an infinite number of edge shapes connecting two vertices that will 
give the same surface area vector as a straight line. These shapes can be called 
equivalent straight line shapes. Thus Eq. (55) is valid for any triangulEU" face with 
equivalent straight line connections. 
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Let a multiply subscripted position vector denote the vectorial average of the 
individual vectors. Thus 

ri23 = g(ri +T 2 +T3) (56) 

is the center of the face. One can easily show that for a plane face with straight 
line edges 

j rdS = Si23Ti2S, (57) 

123 

i.e., the centroid is located at the center. Note that this is not valid if the edges 
deviate from straight lines. From Eq. ( 8 ) it follows that the moment of the area 

Mi 23 = / r X udS = S 123 T 123 X ni23 (58) 

123 

for any face with straight line edges. For a tetrahedral cell with vertices ri, T 2 , rs, 
and T 4 , defined by plane faces and straight edges, the volume is given by 

1^1234 = ^(^2 - ri) X (ra - ri) • (r4 - ri) = ^8123 • (r4 - ri). (59) 

The above formulas can be used to make geometric calculations for an arbitrary 
cell with straight line edges. Each polygonal face can be subdivided into plane 
triangular facets, and the total volume treated as a sum of tetrahedra. The resulting 
surface area vectors and their moments for each face are unique, but the total volume 
will depend on the method of subdivision. Calculations for a regular hexahedral 
cell are presented below. 

A hexahedron defined by eight vertices is shown in Fig. 1, with edges 14, 12, 
arid 15 directed in the positive q, and f directions, respectively. The surface area 
vectors in the positive ^ direction are S 1562 and 84373 , as indicated. The expression 
for the former can be written with the aid of Eq. (55) as 

Si 662 = 2 ^ (60a) 

= (rs6 — ri2) X (ri5 — r2e)- (60b) 

The more efficient form (60a) is expressed as the vector product of the two diagonals, 
showing that each diagonal is perpendicular to the surface normal. The form (60b) 
is in terms of the two vectors joining opposite edge midpoints. Since these vectors 
intersect in the face center ri 562 , it follows that the four edge midpoints and the face 
center all lie in a plane midway between the planes containing the two diagonals. 
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The moment of the area is obtained by dividing the face along one diagonal, and 
can be written as 

Mi562 = 

1562 

In order to insure that the volumes of the hexahedrons sum to the total volume, 
the shape of each face must be precisely defined and consistently used by the two 
neighboring cells. A simple way to calculate the cell volume is to choose an arbitrary 
point inside the cell. The volume is just the sum of the volumes of six pyramids, 
each with one face as the base, and the arbitrary point as the common apex. For 
each non-planar face there exists the location of an equivalent plane face which 
gives the same volume for the pyramid. In fact, there are an infinite number of face 
shapes corresponding to a given equivalent plane face. The volume of the pyramid 
is then one-third of the dot product of the surface area vector of the face and a 
vector from the apex to any point lying in the equivalent plane. 

The earliest expression for the volume of a hexahedron, based on an equivalent 
plane face containing a diagonal, was given by Rizzi (Ref. 9). Unfortunately the 
eqivalent planes for a pair of opposite faces contain oppositely oriented diagonals, 
so that the volumes do not sum properly. Kordulla and Vinokur (Ref. 10) showed 
that of the eight consistent divisions of the faces by diagonals, four result in a very 
simple expression for the volume. If one vertex of a cell main diagonal is chosen 
as the common apex, and the other vertex as the intersection of three equivalent 
plane faces, the six pyramids reduce to three pyramids sharing the main diagonal 
as a common edge. Using the second form of Eq. (59), one obtains the expression 

Vi 2345678 = “(Si485 + S1234 + S1562) ’ (fy ~ Ti). (62) 

Three similar expressions can be derived based on the other three choices for main 
diagonal, each yielding a different value for the volume. 

An alternate expression for the volume is based on the equivalent plane face 
passing through the edge midpoints and face center. Formulas using face centers 
to calculate the volumes of the six pyramids were proposed by Jameson (private 
communication, 1985) and Holmes and Tong (Ref. 11). It was shown by Davies and 
Salmond (Ref. 12) that the same equivalent plane corresponds to a face defined as a 
doubly-ruled surface. Edge midpoints were used to obtain more efficient expressions 
for the cell volume, but full efficiency was not obtained since a cell vertex was chosen 
as the common apex. H the edge midpoint ri 2 is chosen as the common apex, and 
Tqt and r 4 g are each at the intersection of two equivalent plane faces, one obtains 
the expression 

Vi 2345678 = ~[(Ss678 + §2376) * (r67 ~ 1*12) + (S1458 + 84373) ‘ (r48 ~ ri2)]. (63) 


J r X n dS = Ties x 8135 + ri26 x 8126- ( 61 ) 
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It can be shown that the volume given by Eq. (63) is the average of the volumes 
determined by the eight consistent divisions of the cell faces by diagonals. 

For completeness we give the geometric expressions for two-dimensional and ax- 
isymmetric flow. The volume element is now deflned by the co-planar vertices ri, 
T2, Ts, and T4, where k is the unit normal to the plane. For two-dimensional flow 
the surface area vectors take the form 

S12 = (r2 -ri) X k, S23 = (r2 - rs) X k, 

while the area moments become 

Mi2 = j rxndS = ^(ri -ri -r2 -r2)k, 

12 

M23 = r X n ^ (rs • rs - T2 • r2)k. 

23 

The volume is given by 


V"i234 = ^(rs - ri) X (P2 - u) • k. (66) 

For axisymmetric flow, let y be the distance from the axis of symmetry. The 
surface area vectors (per unit circumferential angle) are 

S12 = yi2(r2 -ri) X k, S23 = y23(f2 - rs) X k, (67) 

while the lateral surface area vector becomes 

S1234 = ^(rs-ri) X (r2 -r4). (68) 

The volume (per unit circumferential angle) can be obtained from 

1^1234 = ^[yi32(r3 - ri) X (f 2 - Ti) + yi34(r4 - ri) X (rs - ri)] • k. (69) 
Finite- Volume Discretization. 


Finite- volume methods are based on the discretization of Eq. (1), with the surface 
integral replaced by the sum of integrals over the faces of the cell. The method is 
normally applied to cells defined by the primary grid, so that certain cell faces will 
coincide with the flow region boundary. One can also apply the method to secondary 


( 66 ) 
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cells, in which case the boundary cells are not full cells. The discussion is limited to 
ordered grids, but much of it can be extended to general grids. A full set of integer 
subscripts refers to a cell, or the center or centroid of a cell. Fractional subscripts 
then indicate cell faces, edges, or vertices. Thus the cell is defined by specifying the 
vertices (e.g. these are used to calculate surface area vectors 

(e.g. Sf 1 . . ), area moments, and the cell volume Vf ,• ^ using the formulas derived 

above. This insures that all the geometric identities and constraints are precisely 
satisfied. If F is spatially uniform (which is valid for a uniform free stream), the 
numerical calculation of the surface integral for each face should sum to zero (within 
roundoff errors). This is a necessary condition for the preservation of a free stream. 
Temporal discretization is given by the superscripts n and n + 1, which refer to 
times ti and respectively. 

The primary variable Q is normally associated with the cell, or cell center. In 
some algorithms, such as that of Tong (Ref. 13), it is defined at cell vertices. In the 
steady-state algorithm of Refs. 14 and 15 for the Euler equations, Q is the average 
value on a cell face. Only the first case is treated here. We first consider a flow 
governed by the full set of conservation laws, for which 17 is a function of Q. All 
spatial integrals are replaced by the product of the spatial quantity and the average 
value of the integrand. Thus the geometry of the discretization is treated separately 
from the treatment of the physical variables. There is a certain ambiguity in the 
interpretation of in the relation 

/ U dV ^ (70) 

Strict equality implies that is the average value of 17 in the cell. But in 

order to calculate a surface flux it is convenient to think of as the value of 

U at some average point in the cell, with the = sign replaced by the « sign. This 
interpretation of is also implied when one uses the equation of state to express 
the pressure in terms of the conservative variables. A characteristic of the finite- 
volume method is that the precise location of this average point is not required 
during the calculation. Only in the output of the solution is a location of this point 
desired. Some investigators have suggested the cell centroid for this point. This is 
strictly valid only if all the components of U vary linearly throughout the cell. Since 
the distribution of U is not known, the centroid has no particular advantage over 
the cell center, defined as the vectorial average of the cell vertices. The latter point 
is easier to calculate, and is therefore preferable. Of course, when the finite-volume 
method is applied to a secondary cell, the primary grid point is available as the 
average point for the cell. 

The major problem in a finite-volume method is the calculation of the time in- 
tegral of the average flux over each face. The time integral can be approximated 
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as 


where 






(71) 

(72) 


and At = t 2 — t\. For implicit time integration (9 ^ 0), • t known 

and is therefore linearized about time level n. On the other hand, for a cell face 
that varies with time (Sf . . . . is known and need not be linearized. This is 
one distinction between the finite-difference and finite-volume approach which is 
discussed further in the next section. The calculation of the average flux is now 
presented for each of the class of equations, with comparisons between the finite- 
difference and finite-volume methods made in each case. 


Buler Equations 

For the Euler equations, one only needs the inviscid part of the average flux. 
The numerical procedures for this case can be divided into three classes. In the 
first class y depends only on and but requires an additional 

artificial smoothing flux to stabilize the calculation. This class includes multi-step 
algorithms in which F,.,.! y is alternately F,j,jt and F,+ij,fe, where 

( 73 ) 

and single step centered approximations of the form 

” 2^^*^** (74) 

Another form for the centered approximation is 

= F(77i+ij,*). (75) 

where 

f7.+i,,-,t = i(£7ij,* + t7(+.j',*). (76) 

Form (75) is probably more consistent with the spirit of a finite-volume formulation. 
As discussed in Ref. 11, it gives better results near solid wall boundaries. It has the 
disadvantage of requiring three times as many flux evaluations in three dimensions 
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as does (74). As an alternative, Refs. 11 and 16 suggest a form based on Eq. (2). 
(We are treating the case v = 0 for now.) It can be written as 




(77) 


where ^ 

+ u»+i,y,fe). (78) 

This produces the favorable behavior of form (75) with less computing effort. 

The approximations to the flux integrals for the two faces in the ^ direction based 
on Eq. (74) take the form 



»+^>y.fc 


• + Ft+l,y,fc) 


S- i - fe 


• (F«-i,y,fc + F,-,y,k)]. 


(79) 


The analogous finite-difference expression is 



»+i.y.fc 


Pv+W,*-St 


IJi* 




(80) 


Note that Eq. (80) is simpler than Eq. (79), since it contains fewer terms. Yet 
the conservation that it implies is carried out over a wider region. Specifically, the 
numerical telescoping property is valid over a distance of two cell widths in each 
direction. Thus conservation in a finite-difference discretization is effectively carried 
out over eight sets of overlapping cells with double the size in each direction of the 
original cells. This could lead to larger errors for grids that are not smooth, and 
near flow region boundaries. 

Another source of error in the standard finite-difference discretization for three 
dimensions is the central difference approximation to the surface area vector such 
as y This will result in oscillations for a uniform flow since Eq. (6) is not 
satisfied for the doubly-sized cell. An appropriate difference formula, equivalent to 
the application of Eq. (60b) to the larger cell, was presented in Ref. 17. A more 
efficient formula was given in Ref. 3. Actually, the most efficient form is derived 
from Eq. (60a) as 


S|+i,y,fc - g(rt+i,y+i,A:+i - r,+i,y_i,fc_i) x {vi+ij-i^k+i - r<+i,y+i,fc_i). (81) 

This takes no more operations than the central difference formula and eliminates 
errors for a uniform flow. For problems in which the free stream is a uniform flow, an 
alternate procedure suggested in Ref. 17 is to subtract the free-stream fluxes from 


20 



the conservation equations. This will guarantee exact cancellation of free-stream 
errors resulting from the central difference approximation to y 

To circumvent the need for artificial smoothing fluxes with the first class of ap- 
proximations, upwind-biased approximations that model the waves crossing the face 
are required. For the face with surface area vector Sf , . . = 5^, . , , . . , we 

introduce the notation = nf , , . . • F for the normal flux component. In order 

* T 2 

to achieve higher than first-order spatial accuracy, an upwind-biased numerical flux 
j k states other than and There are two ways this 

is generally accomplished. Probably closer to the finite-volume philosophy is the 
class of approximations associated with the name MUSCL, in which the final cal- 
culation is only in terms of quantities defined at the face. If U~ , . and Ut, . . 

are conservative variables just on the negative and positive sides of the face, then 
the numerical flux is given by an expression of the form 




For first-order accuracy, U. . . = Uij^k- For higher-order spatial accuracy, if 

dimensional splitting is employed, U7 , . . is obtained from upwind-biased inter- 

polation of Ui-ij^ki Uij^k, and possibly I7i+i,y,fc, where the coefficients may be 
modified using appropriate limiters to prevent numerical oscillations. Analogous 
formulas give i . . . Examples of a higher order calculation without dimensional 

splitting for two dimensions may be found in Refs. 18 and 19. 

The solution (82) that represents the wave processes most physically is the one- 
dimensional Riemann solver for the two constant states U7, . . . and U"^, , ... This 
gives a unique value for at the face, and the resultant numerical flux 




(83) 


Since an exact Riemann solver is iterative, and computationally expensive, approxi- 
mate Riemann solvers have been devised, such as those of Godunov (Ref. 20), Osher 
and Solomon (Ref. 21), Collela (Ref. 22), Montagne (Ref. 23), Pandolfi (Ref. 24), 
and Dukowicz (Ref. 25). 

An algebraic approach to obtain (82) is that of flux-vector splitting. For the 
Steger- Warming and Van Leer splittings, the solution is 


k 




(84) 


Calculations using both these splittings for an implicit factored algorithm are pre- 
sented in Ref. 26. Similar calculations using the splitting of Eq. (45) are found 
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in Ref. 6 for an implicit algorithm, and in Refs. 27 and 28 for a two-step explicit 
algorithm. Implicit finite-difference calculations using the Steger- Warming splitting 
are carried out in Ref. 29. Here free-stream fluxes are again subtracted in order to 
eliminate errors in a uniform flow. 


A third approach is based on a local lineairization due to Huang (Ref. 30) and 
Roe (Ref. 5). In terms of the operator some average 

j the solution can be written in the form 






)]. (85) 


The last term can be interpreted as a dissipative correction to a centered approxi- 
mation, with j acting as a numerical viscosity operator. Since (85) violates 
an entropy inequality, the actual form that is used is modified in some manner in or- 
der to prohibit non-physical solutions. j should properly be the Roe average 

of Ut .1 . . and U~, . . , but equally satisfactory resuts are usually obtained if the 
simpler arithmetic average is used. The results of the other two approaches can also 
be writen in a form analogous to (85). For example, Eq. (84) for Steger- Warming 
flux splitting can be written as (85), with the dissipation term replaced by 






( 86 ) 


The other class of higher-order upwind-biased approximations to the numerical 
flux at a cell face involves the consideration of the wave processes at neighboring 
faces. This can be done using any of the three approaches described above to 
represent the wave processes. There is some question as to how the geometry 
of the neighboring face should enter the calculation. From a strict finite-volume 
viewpoint, all calculations relating to the determination of , . . should only 

involve nf , . . , even at neighboring faces. This is computationally expensive, and 

does not tahe into account the grid curvature. It is therefore more appropriate to 
use the local face normal in calculations involving parameters defined at neighboring 
faces. The local geometric scale of the neighboring face can also be involved. To 
illustrate this we consider the application of the local linearization approach. In the 
scheme of Ref. 31, one requires flux differences across neighboring faces, modified 
by appropriate limiters. The relevant parameters defined for face i+ are the 

components of 

(87) 

In the modified flux approach of Harten (Refs. 32 and 33), as implemented for 
curvilinear coordinates in Ref. 34, the relevant parameters are the components of 




( 88 ) 
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where V’j+i is some average of the volumes of the two adjoining cells. The two 
different geometric scales arise naturally from the extension of a Cartesian analysis 
to curvilinear coordinates using Eq. (20). In the finite-volume upwind scheme of 
Coakley (Ref. 35), the parameters are those of Eq. (88) without the geometric 
scale (i.e., - 1.) Further numerical experiments should be conducted to 

determine which geometric scale (or none at all) gives the best results. 

Navier-Stokes Equations 

The calculation of transport terms in the average flux for the Navier-Stokes equa- 
tions necessitates the evaluation of VQ. This can be done using the nonconservative 
expression (17), or basing it on the conservative definition 

j VQdV = j nQdS. (89) 

V s 

The second form leads to a more complex expression, but is more consistent with 
the finite- volume philosophy. It also has some computational advantages, which will 
be indicated presently. Applying Eq. (89) to an auxiliary cell centered at t-|- \ ,j^k, 
one obtains 

^»’+ 2 ) i ,/ c + 2 ^’+2 > Jl ^+ 2 ^* + 2 2 ^’"*"2 2 2 > J >^‘ 

The geometric terms in Eq. (90) can be obtained in two ways. The vertices of the 
auxiliary cell can be defined as vectorial averages of the vertices of the original cells, 
and the exact expressions for areas and volumes can then be applied to them. Since 
true conservation for the auxiliary cells is not required, the averages of the areas 
and volumes of the original cells can be used directly to deflne the corresponding 
quantities for the auxiliary cell. The second method is more efficient, particularly 
for three dimensions. 

The values of Q in the first two terms, corresponding to the longitudinal compo- 
nent of the gradient, are already given. The remaining values, which contribute to 
the transverse components, must be determined by suitable averages. For example 
one can define y+i as 

+ ~ + Q»+l,y+l,*: + Qi,j,k + Qi+l,j,k)‘ (91) 

The sum of the contributions of the last four terms in Eq. (90) using Eq. (91) for 
Cartesian coordinates does not involve Qi,j,k or Qt+i,y,fc- If relaxation methods are 
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used to advance the solution in time implicitly, these terms will affect the diagonal 
dominance of the iteration matrix only weakly due to the nonuniformity of the grid. 
On the other hand if defined as 




(92a) 


or 


^‘+5.y+|.fc " ^(<5t+l,j + l,A: + Qi,j,k), 


(92b) 


the last four terms in Eq. (90) can be combined to improve the diagonal dominance. 
This was pointed out in Ref. 36, where criteria are given on the use of Eq. (92a) or 
(92b). 

The nonconservative form of (V(5),-|.i based on Eq. (17) is usually applied 
in the thin-layer approximation, with the transverse components neglected. Using 
Eq. (18), it can be written as 


(^Q) 








(93) 


with a simple average defining If the transverse terms are required, they 

can be written in several different ways, using expressions such as 




(94) 


with a simple average now defining Note that the transverse terms do not 

affect the diagonal dominance in the nonconservative form. All the above relations 
for the gradient can be used for the transpose by reversing the order of the terms 
in all products. The divergence is obtained by employing a dot product for all the 
products. 

The nonconservative form of a transport term can be related to the analogous 
finite-difference expression. Consider a flux of the form otVQ, where a is a scalar 
function of Q. For the longitudinal component of the fiux integral over face 1 + 5 , j\ k, 
the finite-volume expression is 




(95) 


where (o:/U),.).i defined by a simple average. The analogous finite-difference 
expression is 


( 96 ) 
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It is interesting to compare these two equations (and the corresponding ones for face 
i — \,j,k) with Eqs. (79) and (80) that relate the inviscid flux integrals. We first 
note that the finite- volume expression (95) is now a little simpler. A more significant 
factor is the presence of the volumes in the flux integrals required to calculate the 
gradients. If a central difference approximation is used to determine the volumes 
using Eq. (15), the contributions of the finite-difference flux integrals for cell itj,k 
involve the position vectors r,+2,y,Jt and r»_ 2 ,y,fc. Thus the dependence on the grid 
geometry due to the transport terms is less compact. This could lead to additional 
errors for grids that are not smooth, and near flow region boundaries. One also 
notes that the numerical telescoping of the finite-difference transport flux terms is 
with adjacent cells. It is thus inconsistent with the telescoping of the inviscid terms. 
Since the transport terms give no contribution for a uniform flow, they play no role 
in satisfying Eq. (6). One can thus use Eq. (81) to calculate surface area vectors 
appearing in Eq. (96). Similar conclusions can be obtained for other forms of the 
flux representing transport terms. The transverse component of the flux integral 
has the same behavior as the inviscid flux integral and has the same telescoping 
property. It does exhibit the lack of compactness of the longitudinal component 
due to the presence of the volumes. 

Potential Flow 

Another important case where gradients must be calculated is potential flow. For 
an irrotational flow the velocity is given by 


u = V<6, (97) 

where <f> is the velocity potential. If one further assumes the flow is isentropic, the 
momentum equation has the general Bernoulli integral 

h = C{T)-<i>r--V<l>-V4>^V<j>-rr. (98) 

A 

Here h = e -I- p/p is the specific enthalpy and C{t) is an arbitrary function of t. In 
most problems C is a constant. An arbitrary coordinate transformation has been 
included. The integral is valid for an arbitrary equation of state. For the given 
entropy, the density is some known function 

/> = m- ( 99 ) 

For the special case of a perfect gas f{h) is simply a power of h. Since the state 
of the fluid is completely determined in terms of <f> by Eqs. (97), (98), and (99), <f> 
serves as the single primary variable Q that defines the flow. The conservation of 
mzuss, which has not been used so far, then serves as the single conservation law 
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that determines 4>. We thus have U = p, P = 0 in Eqs. (1) or (13), with G = 0 in 
Eq. (2). Both F and U are now functions of V^. 

There is a fundamental difference in the use of the potential as a primary variable, 
since only its gradient has physical significance. This can be seen by considering a 
very simple solution of the potential equation, namely uniform flow with velocity 
Uoo- The corresponding potential solution is 


<i> = Uoo • r. (100) 

One cannot specify such a flow numerically without defining precisely the locations 
of the position vectors r at which the potential 4> is discretized. While this is an 
obvious statement from a finite-difference viewpoint, in the finite-volume methods 
considered up to now the variables associated with a cell were not precisely localized. 
There is an important consequence in the way gradients are calculated. Applying 
Eq. (17) to Eq. (100) one obtains 


V<A = Uoo • (Ver^ + Vr,T„ + Yfr^). (101) 

The terms inside the parentheses represent the identity tensor analytically. Numer- 
ically, terms such as are approximated the way <f>^ is treated, while the gradients 
of the coordinates are obtained from Eqs. (14), (15), and (18) by approximating 
r^, etc. in some manner. In order to obtain the identity tensor numerically, and 
produce the uniform velocity, the two difference approximations must be the same. 
This result was first obtained in Ref. 37. It also follows from Eq. (98) that for 
a moving grid must be differenced the same way that is treated. In some 
algorithms, the velocity and density may be calculated at different points, therefore 
employing different difference approximations for the derivatives of <f>. The corre- 
sponding derivatives of r used to calculate the gradients of the coordinates must 
therefore also be different. 

From the above discussion it follows that for a finite-volume method, in order 
to preserve a uniform flow , <f> must be defined at precise locations, and the same 
approximations must be used for the derivatives of <f> and the corresponding r in 
calculating V4>. In this respect the finite- volume method for potential flow has 
some of the trappings of a finite-difference method. A clear distinction can still be 
made in the way the cell geometry is handled. If the potentials are defined at cell 
centers, the cell vertices must be defined and used to calculate areas and volumes 
in a finite-volume method. In a finite-difference method, the position vectors at the 
cell centers are differenced in the same manner to calculate both the gradients and 
the geometric quantities appearing in Eq. (20). 

For steady-state algorithms, the areas of the cell faces are the only relevant ge- 
ometric quantities. Two such algorithms are the finite-volume method of Jameson 
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and Caughey (Refs. 38 and 39) , and the implicit algorithm of Holst (Ref. 40), 
as extended by Flores et al (Ref. 37), and Thomas and Holst (Refs. 41 and 42). 
The second algorithm is the basis of the two-dimensional code TAIR and the three- 
dimensional code TWING, both considered finite-difference codes. But according 
to the definitions given above, the finite-volume and finite-difference labels should 
actually be reversed. 

In both methods the conservation law is applied to secondary cells whose vertices 
are obtained as vectorial averages of the primary grid points at which <f> is defined. In 
the method of Jameson and Caughey, all derivatives such as <f>^ and are obtained 
at the secondary grid points by centered box differencing, and are then used to 
calculate fluxes such as at these points. A numerical flux such as . . in 

Eq. (54) is then obtained as the average of the values of F^ at those secondary grid 
points that are the vertices of the i + ^,j,k face. These fluxes are further modified 
by adding recoupling terms to undo the effects of odd-even decoupling, and explicit 
artificial viscosity terms to stabilize the calculation in supersonic regions. The 
above procedure constitutes a finite-difference method according to our definition. 
As is the case for the Euler equations, the area condition (6) is satisfied for two- 
dimensional flow, but is violated for three-dimensional flow. In the latter case the 
box differencing applied to r is equivalent to a secondary grid with piecewise straight 
connections, and each face is the sum of four piecewise facets. The averaging of the 
fluxes at the four vertices does not correctly sum the areas of the four facets for a 
uniform flow. As discussed in Ref. (39), numerical errors for a uniform free stream 
in three dimensions are removed by subtracting free-stream fluxes from the basic 
equation. 

In the TWING code derivatives of <f> and r are obtained at the centers of the 
secondary faces i + etc. by centered differences of known values at the pri- 

mary grid points. These are then used to calculate u and p at the face centers. 
The surface area vectors are calculated from Eq. (60b) in terms of the known posi- 
tions of the vertices. The secondary grid is therefore assumed to have straight line 
connections. The flux integral is then calculated from Eq. (72). Stabilization in 
supersonic regions is accomplished by upwinding the density. In order to minimize 
storage, p is only calculated on $ coordinate faces, where the $ coordinate is chosen 
to be approximately aligned with expected shock surfaces. Its value at other faces is 
obtained by centered averages. The procedure constitutes a finite- volume method 
applied to secondary cells. Note that in two dimensions, the expression for the 
surface area vector as calculated from Eq. (64) in terms of the vertices is identical 
to the corresponding derivative of r at the face center by centered differences. In 
this case the distinction between the two labels disappears. This is due to the fact 
that for steady-state potential flow the flux is only a function of V^, and compact 
differencing is obtained by calculating the fluxes directly on the cell faces. 
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The above discussion maJces it clear that away from boundaries there is a blur- 
ring of distinction between finite-difference and finite-volume methods. In a finite- 
difference method, am auxiliary finite volume grid can be constructed, and the rele- 
vant geometric terms can be calculated from it in order to eliminate those sources 
of error. Conversely, in a finite-volume method it may be necessary to construct an 
auxiliary finite-difference grid in order to calculate gradients properly. 

MOVING Grids 

For an unsteady flow, a grid motion can in general influence the solution of 
conservation laws in three different ways. It affects the convective part of the 
flux due to the presence of the grid velocity in Eq. (2). This is particularly true for 
upwind-biased approximations. The surface area vector over which the flux is acting 
will in general change in both magnitude and direction. Finally, if the grid motion 
is not rigid, the volume of the element will undergo change. In potential flow, there 
is also an additional dependence of p on the grid velocity, as shown by Eqs. (98) 
and (99). In the previous section we saw that in a finite-volume method the grid 
points are treated in a different manner to calculate geometric quantities than the 
way they are used to calculate gradients. Similarly, the temporal treatment of the 
grid points will be different in calculating the change in geometric quantities than 
in the calculation of grid velocities in Eqs. (2) and (98). 

Grid motion relative to a fixed reference frame has been previously studied by 
Thomas and Lombard (Ref. 3). The effect of a non-inertial reference frame has 
been treated by Holmes and Tong (Ref. 11) for constant rotation and an explicit 
integration scheme. A generalization and unification of these results is presented 
here. If a rotating non-inertial reference frame is utilized, one has the choice of 
performing calculations with the fluid velocity components deflned with respect to 
the inertial or non-inertial frame. In the latter instance it is generally assumed that 
external source terms must be present due to the rotation of the non-inertial frame. 
We discuss both situations, and demonstrate that one can perform the numerical 
integration without source terms, thus preserving the strong conservation-law form. 
For concreteness, the discretization of general two-level implicit schemes will be 
presented. 

In the finite-volume method, the grid velocity is treated as a geometric quantity, 
and is interpreted as the rate of displacement of a cell face. In the standard finite- 
difference method, the grid velocity is normally treated as a flow variable, and is 
combined with the fluid velocity to define a transformed velocity. This gives rise to 
unavoidable errors in a uniform free stream. Therefore it is necessary to treat the 
grid velocity in a finite-volume manner to achieve a proper finite-difference method. 
For this reason, the finite-volume formulation will be given in some detail, and 
the relevant changes for the finite-difference formulation will be indicated. Note 
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that the procedures developed in this section are also applicable to space-marching 
algorithms in which the grid changes in the marching direction. 

Formulation of General Grid Motion. 

Let r*(t) be the position vector of a point relative to a non- inertial reference 
frame. Then the corresponding vector relative to an inertial frame has the general 
form 

r(t)=ro(t) + C(t)-r*(0, (102) 

where C is an orthogonal rotation tensor satisfying 

C-C^ = I. (103) 

The absolute velocity v of the point is given by Eq. (3) as the sum of 


Vr = fo + C • r* 


(104) 


and 

Vc = C-f*. (105) 

Here x indicates differentiation of x with respect to time for any quantity x. The 
velocity Vc may be due to the motion of a prescribed boundary surface or the use 
of a flow-adaptive grid. The latter may result from the motion of a free surface or 
the changes in some flow gradients. In the most general situation it could depend 
on all three. 


According to the finite-volume philosophy, the geometric effects due to the grid 
motion are treated separately from the changes in the physical variables. Thus for 
a cell face with surface area vector S(t) we define 


6V = 



n dS dt 


(106) 


to be the volume swept out by the faice during the time interval At. dV^ and dV^ 
are similarly defined in terms of and Vg. Only the sum of the 8Vc for all the 
cell faces contribute to the change in cell volume. The time averaged surface area 
vector for the cell face is defined as 


S = 5n 


if/- 

s(t) 


dSdt, 


(107) 


and the time averaged normal component of the velocity of the cell face is then 


Vn = 6VI{S At). 


(108) 
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It will be convenient to express certain absolute directed vector quantities in the 
non-inertial frame. We thus define 

u* = • u, V* = • V, n* = • n. (109) 

Note that u* and v* are not velocity vectors relative to the non-inertial frame. Let 

B* = • C (110) 

be the antisymmetric angular velocity tensor whose components form the corre- 
sponding angular velocity vector fl*, where 

B*-r* = n*xr*. (Ill) 

If we further let 

Vo = • fo, (112) 

then Eq. (104) can be rewritten in the non-inertial frame as 

v;=vS + n*xr*. (113) 


Discretization with Velocity Expressed in Inertial Frame. 

Assume that the fluid state is given at time n. In general, the position vectors 
r* of the grid points, as well as the quantities ro and C defining the orientation of 
a non-inertial reference frame, will be assumed known at times n and n + 1. Thus 
the grid is assumed to be updated explicitly, even if the flow variables are updated 
implicitly. Unless stated otherwise, we adopt the notation 

Ax = - z” (114) 

and 

x= i(x'* + x'‘+') (115) 

for any quantity x. The spatial indices t,y,A: are omitted for brevity in this section. 

In applying the finite volume method to Eq. (l), the difference in volume integrals 
on the left-hand side can be written as 

A{UV) = AU + AV. (116) 

In evaluating the time integral of the average flux over a face, all geometrically 
defined quantities are assumed held constant at their time-averaged values. The 
flux term can then be approximated as 

[(1 - ^)F” + At, (117) 
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where the quantities ii and Vn defined by Eqs. (107) and (108) are used in evaluating 
Fn = n • F at both time levels. These same quantities are also used in evaluating 
the flux Jacobian A matrix when time linearization is used. 

It is simpler to calculate the geometric quantities in the non-inertial reference 
frame. From the known position vectors of the cell vertices, the surface area vectors 
S*, area moments M*, and cell volumes V can be obtained from Eqs. (60) through 
(63) at times n and n + 1, and their time averaged values calculated from Eq. (115). 
(F is required if a production term is present in Eq. (1).) The conservation of 
volume condition (9) can be satisfied by calculating 6Vc for each face exactly from 
the known positions of the four vertices at times n and n + 1. A simpler method 
is to obtain an approximate displacement 6t* as the average of the Ar* of the four 
vertices, and calculating SVc from 


(118) 

One must then use the sum of the SVc of file six faces in place of the true AV in Eq. 
(116). The simpler method should normally suffice, although second-order accurate 
algorithms with very large grid distortions may require the exact procedure. 

If the grid is given in a non-inertial reference frame, one obtains SVr for each face 
from 

SVr = (vS • s* + n* • M*) At. (119) 

Since the S* and M* satisfy Eqs. (6) and (8) when summed over all the faces, the 
SVr will sum to zero, as required by rigid body rotation. The angular velocity n* 
is obtained from the components of B* calculated as 

B* = • AC. (120) 

One can readily show that Eq. (120) results in an antisymmetric tensor. The trans- 
lational velocity Vq is similarly calculated as 

vS = • Aro. (121) 

The surface area vector S in the inertial frame can be calculated from 

S = -[C'*-(S*)'‘ + C'‘+^ •(S*)'‘+^] (122) 

2 

or 

S = C • S*. (123) 
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The second expression is simpler, and differs from the first by an error of order 
(At)^. Note that for both expression S will differ from 5* by that same order of 
error. 

The final step is the calculation of from Eq. (108), where 

6V = 6Vc + 6Vr. (124) 

Any finite-volume algorithm discussed in the previous section can now be applied, 
with the grid motion included in the calculation of the flux. The change in cell 
volume due to grid distortion manifests itself as an additional explicit term, as 
shown by Eq. (116). 

In a finite-difference method, the positions of the cell centers are known at times 
n and n + 1. Since the grid motion just affects the convection part of the flux, it 
is sufficient to examine the inviscid flux terms only. For central-difference approxi- 
mations, the integration of Eq. (20) over a time interval results in expressions such 
as Eq. (80). The evaluation of geometric quantities at each grid point is done in 
a manner analogous to that in the finite volume method, taking into account that 
conservation is effectively applied to a doubly sized cell. Condition (21) is satisfied 
by calculating the surface area vectors from Eq. (81) in terms of the positions of the 
vertices of the doubly sized face passing through the grid point. The area moment 
can be analogously evaluated by applying Eq. (61) to those vertices. 

Since the central-difference approximations for the effective cell volumes do not 
define a precise volume, one must obtain an approximate 6Vc from Eq. (118), with 
5r* replaced simply by Ar* at that grid point. (It is possible to calculate a volume 
for the doubly sized cell that is consistent with the area formula (81), using either 
Eq. (62) or (63). With this additional complexity, any advantages of the finite- 
difference discretization are lost, and it is better to go directly to the more compact 
finite-volume discretization.) The summation of the 6Vc (with appropriate signs) 
for the six faces of the doubly sized cell to obtain the appropriate AV to use in 
Eq. (116) is equivalent to solving Eq. (22) by central differences. The importance 
of using Eq. (22) was first pointed out in Ref. 3. The effects due to a non-inertial 
reference frame are handled precisely as in the finite-volume case. 

Discretization with Velocity Expressed in Non-Inertial Frame. 

For many applications in which the grid is defined in a non-inertial reference 
frame, it is more convenient to employ velocity components referred to axes fixed 
in that frame. The momentum conservation law for these components requires the 
presence of source terms. There is an analogy with the use of curvilinear coordinates, 
where the momentum equations for the curvilinear velocity components introduce 
source terms. But these terms can be eliminated by writing the equations for the 
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Cartesian velocity components instead. Similarly, by employing components of the 
absolute velocity (as expressed in the non- inertial frame), one can also eliminate 
the source terms. 

Let the column vectors U* and F*, and the matrix A* be defined in terms of u*, 
V* and n* given by Eq. (109). Using the rotation matrices 



1 

0 

0* 


1 

0 

o' 

c = 

0 

c 

0 

, c-i = 

0 

0 


0 

0 

1 


0 

0 

1 


one can write the transformations 

u = cu‘, f„ = cf;i, a = ca‘c-K 


(125) 


(128) 


In order to reproduce a uniform flow in the inertial frame using an impiicit scheme, 
one must use an incremented quantity proportional to AU, where 


This suggests defining 
and 


AU = C^+^AU* + ACU*^. 


AU* = {C~^)^+^AU 


B* = ■^(C-^)”+^AC = 


0 0 0 
0 B* 0 
0 0 0 


where 


B’ = ^(C’’)’>+‘ - AC. 
At' ' 


(127) 


(128) 


(129) 


(130) 


Note that B* differs from B* by a term of order (At), and is therefore not anti- 
symmetric. Equation (127) can now be rewritten as 


AU* = AU* - B*U*^ At. 


( 131 ) 


The second term is only present in the momentum equation, where it represents a 
numerical Coriolis correction. 

Before writing the implicit equation for AU* , we define the tensor 

i= (C^)”+^ -C” =I-B* At, (132) 
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and the corresponding matrix operator 






10 0 
0 10 . 
0 0 1 


(133) 


Note that they differ from the identity tensor and unit matrix, respectively, by terms 
of order (At). Premultiplying Eq. (116) by (C"^)”"^^, we obtain for the difference 
in volume integrals the expression 


V^+^AU*+iU*^ AV. 


(134) 


The corresponding flux integral can be written in the time linearized form as 


{IF* + 6A*^AU*)S* At. (135) 

Equations (134) and (135) define an implicit algorithm to calculate AU*. This is 
then corrected by the Coriolis term using Eq. (131). The scheme is fully conserva- 
tive, and preserves a uniform free stream. Note that in order to achieve the strong 
conservation-law form it is necessary to use the dot product of the explicit term in 
the momentum equation with I, and to use the proper numerical representation of 
the Coriolis term. 

The fluid velocity vector relative to the non-inertial frame is given by 

u*' = u*-v*. (136) 

In order to calculate it one requires v* at cell centers at times n and n 4- 1. The 
positions of the cell centers are given for a finite-difference method. For the finite- 
volume method one must define the location of the cell center when the non-inertial 
frame undergoes rotation. This situation is similar to the one discussed previously 
for potential flow. In addition, in order to achieve second order accuracy in time, 
one must know Tq, C, and r* for each grid point at three time levels. Thus, using 
Eq. (136), one can treat the relative motion in an arbitrary non-inertial reference 
frame in a completely conservative way. 

Treatment of Boundaries 

The treatment of flow region boundaries depends on whether conservation is ap- 
plied to cells defined by the primary or secondary grid. These will be referred to as 
finite-volume and finite- difference grids, respectively, even though the finite-volume 
and finite-difference methods can be applied to grids of the other family. The most 
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physically satisfying way of treating boundaries is based on characteristics consider- 
ations. This is a natural procedure for upwind methods, and the most appropriate 
for open boundaries. Since this procedure is fully described by Chakravarthy (to 
be published), we discuss only those that use the basic equations as applied to wall 
boundary conditions. These can be conveniently divided into two main classes. In 
one class the unknown variables on the wall are integrated together with those at 
interior points by an appropriate application of the conservation laws and bound- 
ary conditions. In the other class, the unknown boundary quantities needed to 
calculate flux terms at interior points are obtained using extrapolation, reflection 
principles, or some auxiliary equation. We will show how the conservation laws 
can also be utilized in this second approach. Both approaches will be discussed for 
flnite-difference and flnite-volume grids. Note that questions of stability and pro- 
gramming efficiency, although both very important, are beyond the scope of this 
paper. The type of grid also plays an important role at zonal boundaries and in the 
treatment of grid singularities. These topics axe also treated in this section. Only 
stationary grids will be considered for simplicity. 


Wall Boundary Conditions for Finite-Difference Grid. 

Let the wall be a constant f surface with index k = 1. Associated with the 
boundary point t,i, 1 is a secondary grid half-cell whose center is designated as 
the point i,j, |, as shown in Fig. 2. Since the points r,-^,o do not exist, the finite 
difference expressions for and V previously derived for interior points are not 
valid at the points i,j, 1. These quantities are required to calculate transport terms 
for points t,y, 2 and t,j, 1, as well as the geometry of the boundary half-cells. The 
simplest way to modify the interior formulas is to replace r,-,y,o by and multiply 
the final result by 2. As an example, formula (81) becomes 



1 

4 


(r<,y+i,2 




(137) 


The use of first-order instead of second-order accurate one-sided expressions is ac- 
tually preferable for two reasons. It prevents possible unphysical answers (e.g. 
negative volumes) for grids with spacing discontinuities. It is also necessary in or- 
der to satisfy Eq. (6) for the doubly-sized half-cell associated with the boundary 
point. When calculating transport terms at point i,j,l using Eqs. (17) and (18), 
one should use 

= iQw - (138) 

for consistency. Therefore the point i,j,2 should be chosen close enough to the 
boundary to be in the linear range of the variation of Q when transport terms are 
important. 

Extending the work of Thomas (Ref. 43), we apply Eqs. (20) and (54) to the 
half-cell centered at the point t,y, |. Second-order accurate spatial differencing 
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yields 


+ (Phui - = «• 


(139) 


where 


Ff 


A /3„ 1„ ,,, 

,3 


(140) 


4®'*+i 


Note that a partial finite-volume approximation is employed. Thus we use Vi-y. 1 and 
y 1 (instead of their values at A; = |), since these have already been defined so 
as to satisfy the geometric conditions. 

We can eliminate f/»,y ,2 using the conservative difference approximation at point 
t,y, 2. Introducing 

P (141) 

we obtain for the second order accurate conservative spatial differencing of Eq. (20) 
at point t,y, 1 the relation 


+ (« + ( 142 ) 

+ - PLu>)/^( + = 0 . 


where 




(143) 


An important special case of Eq. (142) results from the solid wall condition 


(Sf -m),-,,,! =0, (144) 

valid for all flows. By combining Eq. (144) with the dot product of S^ y j and the 
second component of Eq. (142), one obtains a time- independent relation expressing 
the conservation of normal momentum at the boundary point. 

Consider first the solution of the Navier-Stokes equations for a perfect gas. The 
no slip condition 

m,,y.i = 0, (145) 

also results in 

P.,y,i =Ke„y,i. (146) 
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If the wall temperature is prescribed, and is the wall specific internal energy, 
then we also have the condition 


c*,y,i — ( 147 ) 

Thus there is only one unknown variable on the wall. 

For the first class of methods, Thomas (Ref. 43) suggests using the continuity 
equation, since the boundary conditions supplant the momentum and energy con- 
servation laws. The first component of Eq. (139) or (142) would thus be the missing 
equation to use in an implicit algorithm. Due to the presence of in Eq. (139), 
it is perhaps more appropriate to use Eq. (142), since the boundary conditions only 
specify the time derivatives of Actually, the wall temperature condition only 

relates to pi,j,\ through Eq. (147), so that one could integrate the energy equa- 
tion instead. Whichever equation is chosen, conservation for the remaining conser- 
vative variables will not be precisely satisfied for the boundary half-cell. Therefore 
some interior numerical flux terms will not be precisely cancelled, and conservation 
of some variables for the entire flow region will be violated in an integral sense. 

As pointed out by Mehta (private communication, 1985) any independent relation 
can be used to integrate with time. The time-differenced form of the normal 
momentum equation obtained from the second component of Eq. (142) provides such 
a relation for the pressure. This can then be used instead of the energy equation, in 
view of Eq. (146). The use of the time derivative of a conservation law, as opposed 
to its direct employment when solving the continuity or energy equation, can have 
important implications in the implementation of a factored implicit algorithm. This 
is discussed further below. 

For the second class of methods, the pressure is the only unknown wall quantity 
needed to calculate flux terms when applying the conservation laws at interior points 
i,j,2. The time-independent conservative normal momentum equation at the point 
i,y, 1 provides a relationship between the wall pressures and quantities at interior 
points. In implicit algorithms both the direct and time-differenced forms of this 
relationship are required. In this respect the second method has aspects of both 
versions of the first method. 

The practical implementation of all of these approaches could require further 
approximations which decrease the spatial or temporal accurax:y of the algorithm 
at the boundary, and may involve a restriction to orthogonal grids. A factored, 
implicit, central-differenced implementation of Eq. (139) can only be first order 
accurate in time due to the presence of 17*, y, 2 - Flux terms such as y 2 must 
be treated explicitly, thus again reducing the temporal accuracy to first order. One 
can maintain second order accuracy in time by letting — ^»,y,i y 2 = 

F,.,.! y 1 , but then the spatial accuracy is reduced to first order. If one uses Eq. (142) 
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instead of (139), one must still treat some flux terms explicitly and suffer a loss in 
temporal accuracy. 

Another problem arises with the use of the time-differenced normal momentum 
equation for large wall curvature or severe non-orthogonality in the grid. Under 
those circumstances the pressure change at point i,j,l is coupled to the pressure 
changes at neighboring points on the boundary as well as the points t, j, 2 and t, jf, 3. 
In an implicit updating of using this equation, all spatial difference operators 
operating on are of order one. Consequently a factored implicit algorithm 

is not possible. Similarly, in the second method a factored implicit treatment of the 
boundary flux for point i,j,2 cannot be accomplished. Thus the use of the normal 
momentum equation in a boundary procedure for factored implicit algorithms is in 
principle limited to moderate wall curvature and only small non-orthogonality in 
the grid. Fortunately in most practical situations the grid spaxring in the direction 
is much smaller than the grid spacings along the boundary. The resulting coupling 
of pressure changes at points on the surface is negligible, and the above restrictions 
can be removed. 

If the heat flux Qu) is prescribed at the wall instead of the temperature, condition 
(147) is replaced by 

(/jn • V£)j^y^l = Cv9u;> (^^®) 

where k{e) is the thermal conductivity and is the constant speciflc heat at con- 
stant volume. Using Eq. (145) one has the identity 


(V£).j,i = (149) 

For steady-flow calculations. Hung and MacCormack (Ref. 44) suggest using VJST to 
evaluate the wall heat flux, since the total enthalpy H has a more linear behavior 
than £. For unsteady flow, it is probably more appropriate to use V(e/p). In 
the first method of satisfying the boundary conditions, one now integrates both 
the continuity and energy equation at point *,i, 1 using either Eq. (139) or (142). 
The alternative approach uses the time-differenced forms of the normal momentum 
equation and Eq. (148). These two equations are also used in the second method 
to calculate boundary flux terms when updating U»-,y, 2 > 

For the solution of the Euler equations the only boundary condition is Eq. (144). 
In the first method of treating the boundary one would now integrate the continuity, 
energy, and tangential momentum equations at point t,i, 1 using either Eq. (139) 
or (142). Alternatively one could use the time-differenced form of the normal mo- 
mentum equation instead of the continuity or energy equation. Since the pressure 
is the only unknown required to evaluate boundary flux terms, one would only need 
the normal momentum equation to apply the second method. However if the wall 
curvature is not negligible, the expression for the pressure at point i,j, 1 involves the 
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velocity and density (as well as the pressure) at neighboring points on the bound- 
ary. For example, the central difference approximation to the normal momentum 
equation using Eq, (142) produces the boundary terms 




(8-^)(S* -S'p),-,,-,! - •S^,)(pu-S*)j+lj,l +S?j. j •(pS^),+l,,-,l| + - • • 

(150) 


There are several ways to evaluate these new unknown quantities. The simplest is 
to replace them by their values at the closest interior points, with an attendant loss 
of spatial accuracy. A more accurate procedure is extrapolation of interior values 
to the boundary. In the spirit of the conservation laws it is more appropriate to 
extrapolate p and m = pu, with the latter then projected onto the surface to satisfy 
Eq. (144). If the grid is highly non-orthogonal, extrapolation should probably be 
done in the direction normal to the surface, rather than along the f coordinate. A 
third possibility is to utilize Eq. (139) or (142) to update p and m on the surface 
after the interior points have been updated. This is the most rational procedure, and 
should probably be used in the calculation of highly unsteady flows. As discussed 
earlier, an implicit treatment of the pressure at point t,i, 1 is not feasible in a 
factored implicit algorithm for large wall curvature or severe grid non-orthogonality. 


For potential flow there is only one conservation law, and Eq. (144) gives the 
required boundary condition. Due to the manner in which the velocity potential <f> 
on the boundary is used in solving the equation at point one must use the 

first method to satisfy Eq. (144). The proper procedure is to solve Eq. (139). This 
is more accurate than using simple reflection principles. 


Wall Boundary Conditions for Finite- Volume Grid. 

Let the wall be a constant ^ surface and the boundary cell be designated as i,j, 1. 
The boundary face then has the designation t, j, The position vectors 




(151) 


and vertices of a boundary half-cell whose center is designated 

as the point i,j, |, as shown in Fig. 3. In order to calculate transport terms and 
describe the geometry of the boundary half-cell we require the quantities Sjy j, 

i, and ^ Sjy ^ is defined by Eq. (60) in terms of the vectors r,±i 

and y are defined as twice the corresponding expressions for the half- 

cell. Thus 


(r, 


‘.j-i 


(152) 
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When calculating transport terms at points »,y, 1 and *,y, | using Eqs. (17) and 
(18), one should use 


~ ~ 2(Q»,y,i 


(153) 


Applying Eq. (5) to the boundary half-cell centered at the point i,j, |, one obtains 
the second-order accurate expression 


where 


+Fi+ij,i),etc. 


(154) 


(155) 


Note that we use and ^ (instead of their values at ik = |), since these 

have already been defined so as to satisfy the geometric conditions. 

We can eliminate C^»,y,i using the finite volume approximation for cell t,y, 1. In- 
troducing 

^ = Vij.ilVi.i,x, (156) 

we obtain for point t,y, ^ the second order accurate conservative relation 


• F)*,y,a + 4(S' . F).., - (4 - ^)(S= ■ F).„. r 


(157) 


where 


»+2 2 




^)»+§,y,i 




• Fj ^.1 y i,etc. (158) 


The treatment of wall boundary conditions for finite-volume grids basically follows 
that described for finite-difference grids, with Eqs. (139) and (142) replaced by 
Eqs. (154) and (157), and the finite-difference boundary point »,y, 1 replaced by 
the finite- volume boundary point t,y, Note that conservation of all conservative 
variables for the entire flow region is automatically satisfied when a finite-volume 
grid is used. The purpose of introducing a boundary half-cell in this case is to relate 
boundary values to interior values in a conservative manner. 
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One can again apply the first method of treating the boundary. Note that in 
this case there are two boundary cells, both sharing a common boundary, that need 
special treatment. The boundary variable Uij^i now represents an average value of 
U over a cell face instead of a cell-averaged quantity. The points t, j, | and t, j, 1 are 
a half-cell width apart. All these factors make an implicit algorithm more involved 
near the boundary. Nevertheless it is probably worth pursuing if one requires exact 
conservation in the integral sense and accurate representation of all flow variables 
on the boundary. 

The second method of treating the boundary is obviously the natural one for a 
finite-volume grid. The normal momentum equation based on Eq. (157) can be 
utilized to various degrees of approximation for this purpose. Since the point t , j, 1 
is closer to the boundary than the corresponding point for a finite-difference grid, 
one would expect these approximations to be more accurate for the finite-volume 
grid. 

Zonal Boundaries. 

Another situation where the difference in the type of grid is important is the 
case where the boundary is a zonal boundary between two regions with completely 
disparate grids. The two types of grids are illustrated in Fig. 4 for two dimensional 
flow. For the finite-difference grid, the dots show the location of the grid points in 
one zone, and the squares the grid points in the other zone. A conservative zonal 
boundary procedure requires the interpolation of data between the two grids on the 
zonal boundary, and the partitioning of flux on a flux conservation line arbitrarily 
chosen to lie in one of the two zones. Special handling of the fluxes is required on 
sides AFE and BCD of the cell straddling the zonal boundary. Excellent results 
have been obtained in flow calculations using this procedure by Rai (Refs. 45-47). A 
finite- volume grid adapts more naturally to the zonal boundary, as shown in Fig. 4b. 
The partitioning of the flux can now be carried out directly on the zonal boundary, 
leading to a conceptually simpler algorithm. Calculations using a finite-volume 
grid have been carried out by Eriksson and Rai (to be published) and Walters et al 
(Ref. 48). 

For both types of grids, the above mentioned zonal procedures carry out the 
partitioning of the flux in terms of the normal flux component F’r. As pointed out 
in Refs. 45 and 48, if the flux conservation line is curved one cannot simultaneously 
conserve the flux and maintain a uniform free stream. However, for a finite-volume 
grid an alternate procedure is possible which can accomplish both objectives. For 
a given boundary cell one can define a separate boundary face for eax:h cell in the 
other region. Thus the zonal boundary of cell 1 in Fig. 4 consists of the three faces 
FE, ED and DC. If a unique flux is assigned to each boundary face, then both flux 
conservation and free-stream maintenance is possible. For first-order accuracy, the 
flux across FE is calculated from U\ and U 2 ', the flux across ED from U\ and 1 / 3 ; 
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and the flux across DC from U\ and . Higher order flux calculations will involve 
more complex dependence on the values of U in neighboring cells. The MUSCL 
approach, which was also used in Ref. 48, is probably the best one in this situation. 
Note that boundary faces for three-dimensional zonal boundaries will no longer be 
quadrilaterals in general, and expressions for a general polygonal face must be used 
to calculate surface area vectors. For time-accurate calculations in two and three 
dimensions, more general formulas for volumes of polygons and polyhedra must be 
used, since the boundary cells will in general not be quadrilaterals or hexahedra. 

Grid Singularities. 

There are two types of grid singularities that can occur on a boundary. One 
is due to a physical corner occurring on a solid boundary that is described by a 
single coordinate surface. The primary grid point deflning the corner is called a 
real singular point. The other type occurs when a primary grid point on a smooth 
solid boundary deflnes a corner in a coordinate surface. Such a point is termed 
a topological singular point. For a finite-volume grid these singular points are 
cell vertices, and their singular nature will not affect the evaluation of geometric 
quantities based on straight line connections. Since the finite volume algorithms 
do not involve flow variables at the cell vertices, the algorithms will also not be 
affected by real singularities. Some modifications in the boundary procedure may 
be necessary for cells whose vertices include topological singular points. In either 
case a large loss in spatial accuracy can be expected, since the singular nature of 
the grid or the flow is ignored. 

For a finite-difference grid both U and geometric quantities must be defined for 
the singular points. The non-analytic behavior of U at a real singularity cannot be 
simply expressed. Consequently it would be difficult to account for it properly in a 
numerical algorithm. On the other hand, the non-analytic nature of the grid at a 
singular point can be expressed algebraically. This knowledge has already been uti- 
lized to generate algebraic grids with singular corners in Ref. 49. A similar approach 
can be used to modify a finite-difference algorithm near a topological singularity. 
We illustrate the procedure for a Navier-Stokes central-difference algorithm near an 
H-type singularity in two-dimensional flow. 

Figure 5 depicts the grid in the neighborhood of the H-type singularity at ^ = 0 
and r/ = 0. The part of the coordinate line r; = 0 for positive $ lies on the solid 
body, while the part for negative ^ lies in the flow. We assume the two paxts are 
normal to each other at the singularity. The second method of applying boundary 
conditions will be used, with the wall curvature neglected. All calculations that 
involve quantities defined at point 0,0 will have to be modified. The geometric 
quantities defined in Eqs. (14) and (15) reduce to 

= r,, X k. S’' = k X r^, V = • r,, x k (159) 
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for two dimensions, where k is the unit normal to the plane. 

Consider first calculations for the grid point 0, 1. The grid is non-analytic between 
this point and the singularity. One can easily show (see Ref. 49) that to a good 
approximation 

r(0, ri) « ro,o + (ro.i - ro,o) (160) 

and 


s^o,n) 


® 0,1 






’ 0,1 




'"(0,»7) 


Vo.i 

rilin' 


(161) 


between these two points. While at point 0, 1 is calculated by the ordinary 
central difference formula, r,, is calculated as 


(r„)o,i = (ro,| - ro,i)/A»/. (162) 

With Tq I obtained from linear interpolation, and Fq i obtained from Eq. (160), 
one can rewrite Eq. (162) as 

(r„)o,i = - (V^ - l)ro,i - (2 - V^)ro,o]- (163) 

The geometric quantities Sq i, Sq i, Vb,i, (V^)o,i, and (Vr;)o,i are then obtained 
from Eqs. (159) and (18). 

The discretization of Eq. (20) at point 0, 1 requires modification of the term 

iP:i)o,i = (^ 0 % - ( 164 ) 

can be calculated by the standard methods in terms of quantities already 
defined. The variation of the flow variables with r is analytic. Therefore between 
the points 0, 0 and 0, 1 we can expect the inviscid part of the flux to satisfy 


F(0,»;) « Fo,o + (Fo,i — Fo,o)x/»t7a^. (1®^) 

Using Eqs. (161) and (165) to evaluate we can write Eq. (164) as 

for the inviscid part of the flux. The transverse component of the transport part of 
the flux has the same behavior as the inviscid part. In calculating the longitudinal 
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component, we note from Eqs. (18) and (161) that is approximately constant 

between the points 0,0 and 0,1. Assuming that the dependence of both a and Q 
on T] is given by Eq. (165), one obtains for the longitudinal transport flux term at 
0, 1 the relation 

(aS” • VnQn)o,i = 1“°.' + (\^ - l)“o,o|(S'' ■ SVt'jo.iWo.i - «o,o)/Aij. (167) 

The calculations for the grid point —1,0 follow the same procedure as for the 
point 0,1, with the roles of r) and i reversed. At the point 1,0 the only quantity 
one requires that must be modifled is 

(r^)i,o = ^^[^ 2,0 - (V^ - l)ri,o - (2 - \/2)ro,o]. (168) 

The only quantity required at the singular point 0,0 is the normal gradient of Q 
which enters into boundary conditions. It is calculated as 

(n • VQ)o,o = (Q- 1,0 - Qo,o)/|r_i,o - ro,o|- (169) 

Note that geometric quantities at the singular point are undefined, eind are therefore 
never used. If the first method of satisfying boundary conditions is employed, one 
could perform calculations for the half-polygonal cell associated with the singular 
point, using the concepts described above. 

All the above formulas assume that the grid is sufficiently smooth. If the grid is 
generated algebraically, the techniques presented in Ref. 49 can be used to guarantee 
sufficient smoothness. In numerical grid generation schemes based on the solution 
of partial differential equations, the corner in the »; = 0 curve will propagate into the 
interior using the standard procedures. One must therefore use the above concepts 
to modify the numerical grid generation scheme. This involves using approximations 
such as Eqs. (160) and (161) to derive appropriate difference approximations to 
partial derivatives similar to Eqs. (163) and (168). 

Strong and Weak Conservation 

The analysis in this paper has been carried out using vector notation, even though 
all computations are ultimately performed in terms of scalar components. This was 
done for two different reasons. The formulation in terms of physical vectors is more 
compact, and gives more emphjisis to the physical content of the numerical methods. 
Secondly, there are a number of different ways to obtain scalar equations from the 
vector equations. Most of them are motivated by a desire to keep the equations in 
strong conservation-law form. This terminology was coined by this author in his 
earlier paper devoted to the differential formulation of conservation laws (Ref. 1). It 
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refers to the form of Eq. (13) in the absence of true, physical source terms (P = 0), 
in which all terms are derivatives with respect to the independent variables. The 
decomposition of the vector momentum equation into scalar equations in ter ms of 
curvilinear velocity components results in additional undifferentiated terms that act 
like fictitious sources. This has been termed a weak conservation-law form by the 
author. Various methods of obtaining scalar equations in strong conservation-law 
form are discussed in Ref. 1. The one most often used in practice is to write the 
Cartesian components of the vector conservation law. Actually, in some algorithms 
(e.g. Ref. 50) it is advantageous to employ the weak conservation-law form. 

In the section treating moving grids we encountered another situation that is 
normally thought to require the weak conservation-law form — namely the use of 
a non-inertial reference frame. It was shown there how the strong form can be 
preserved by using the absolute velocity. There are two other situations that re- 
sult from ignoring a space dimension in which the equations are usually expressed 
in weak conservation-law form. They are quasi-one-dimensional and axisymmetric 
flow. We now show that in a proper finite-volume formulation, the undifferenti- 
ated terms become boundary terms for the ignored direction. Thus the integral 
conservation law can always be satisfied for these cases. 

Quasi-One-Dimensional Flow. 

The differential formulation of quasi-one-dimensional flow results from apply- 
ing the coordinate transformation x(^) to a one-dimensional channel whose cross- 
sectional area vector is 

= 5(x)i, (170) 

where i is the unit vector in the x-direction. In terms of the cell volume V = Sx( 
and the velocity component u = i-u, the continuity and energy equations for inviscid 
flow become 

(pV)r + (puS)( = 0 (171) 

and 

(eV)r + [(c +p)u5]^ = 0. (172) 

Since there is no flux of mass or energy at the wall, these two equations are in strong 
conservation-law form. The momentum equation is usually written as 

(puV)r + [(p + />u^)<S]^ — = 0. (173) 

The source term pS( results from the pressure acting on the channel walls. 

In order to circumvent the weak conservation-law form, some investigators (see 
Ref. 51) write the momentum equation in the quasi-conservative form 

(puV)r + (pu^5)^-f-5p^ =0, (174) 
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in which differentiated ternas are multiplied by geometric quantities. While this 
form is strictly not conservative, it satisfies the weak form of the jump conditions 
(11), and should possess the appropriate shock capturing capability. On the other 
hand, the weak form is only apparent, not real. A finite-volume discretization of 
the momentum conservation law for cell t gives 

Vi{pu)ri + [(p + pu^)S]i^i - [(p + ptt^)5],_i - p«,t(5i^.i - Si_i) = 0. (175) 

Here pty, is the average pressure acting on the cell wall, while p, is the average value 
of p throughout the cell. Thus the source term is actually a boundary term. The 
boundary condition that relates the two pressures is derived from the transverse 
momentum equation. According to the quasi-one-dimensional approximation, this 
relation is 

Pwi « Pf (176) 

Thus the undifferentiated term only appears to be an interior source term when 
relation (176) is used as an exact identity to eliminate p^i in Eq. (175). When one 
realizes that Eq. (176) is only an approximate relation, resulting from the approxi- 
mate solution of the transverse momentum equation, one sees that the last term in 
Eq. (175) is really a transverse boundary term. It should be treated in an algorithm 
in the same manner as any other boundary term. 

The proper way to discretize the equations is to apply the conservation laws to 
the primary grid cells, which are defined by specifying the values of and 
at the cell boundaries. In a finite-volume discretization, one also needs the cell 
volumes, which axe given by 

Vi = \{xi+L - x,-_i)(5i^.i -h 5,_i). (177) 

Due to the quasi-one-dimensional approximation, the above relation is valid for 
both planar and axisymmetric flow. One can also apply the finite-volume method 
to secondary grid cells, with special treatment of the half-cells at the two ends of 
the channel. 

Axisymmetric Flow. 

The differential formulation of axisymmetric flow differs from the two-dimensional 
case primarily in the definition of geometric quantities. If k is the unit normal to 
the rj plane, and y is the distance from the axis of symmetry, then one uses 


for two-dimensional flow and 


r,=k 

(178) 

II 

(179) 
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for axisymmetric flow in Eqs. (14)and (15). Since = 5^k, one can also write 
V = for two-dimensional flow and V = yS^ for axisymmetric flow. The strong 
conservation-law form (20) with P = 0 and the f derivative term absent results for 
both flows, except for the momentum equation. If j is the unit normal to the axis 
of symmetry in the t) plane, and t; = j • u, then the momentum equation in the j 
direction for axisymmetric, inviscid flow is usually written as 

{pvyS% + (F«) J + - pS< = 0. (180) 


The weak conservation-law form of Eq. (180) is again only apparent. The conser- 
vation law should really be applied to a wedge-shaped region of angular width A(p, 
where <p is the circumferential angle. Each cell extends between the two planar 
faces of the wedge, which act zis cell boundaries in the circumferential direction. 
There is no convection of fluid through these boundaries, and the only contribution 
to the flux is the pressure in the momentum equation. The circumferential com- 
ponent of that equation is identically satisfled by axial symmetry. A finite-volume 
discretization of the radial component for cell i,j gives 
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(181) 


The last term represents the contributions from the two circumferential boundaries, 
where (pb)t,j is the average pressure acting on the lateral faces of cell t, j. Terms such 
as y are defined as in Eq. (72), while the geometric quantities are calculated 
from Eqs. (67)-(69). The condition of axial symmetry yields the boundary condition 

(Pb).-,j » Pi,r (182) 


Note that this relation is not exact, since the average pressure acting on a lateral 
surface does not equal the pressure determined by the conservative variables av- 
eraged over the cell. Since the angular width A<p is arbitrary, we can make the 
further approximation 

sin « 1. (183) 

Actually, by choosing A(p small enough, one can make the error in Eq. (183) less 
than the round-off error in the computation, so that Eq. (183) becomes effectively 
an identity. We thus again find that the undifferentiated term only appears to be 
an interior source term when relations (182) and (183) are used in Eq. (181). 

The proper way to discretize the equations is to apply the conservation laws to 
primary grid cells. The axis of symmetry then serves as a boundary of zero surfaw:e 
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area for a row of wedged-shaped cells, and consequently does not contribute to any 
flux calculations for those cells. One thus avoids any difficulties due to the axis 
singularity. 

In summary, we have demonstrated for both quasi-one-dimensional and axisym- 
metric flow that a finite-volume discretization applied to primary grid cells enables 
us to preserve the strong conservation-law form. The undifferentiated terms are 
actually boundary terms for the flow region, and should be handled in the same 
manner as other flow region boundary terms. While the axisymmetric case was de- 
rived for an inviscid flow, it can be readily extended to a viscous flow by including 
the normal viscous stress in calculating the force acting on the lateral surface. 

Concluding Remarks 

This survey of finite-difference and finite-volume approaches has revealed that 
comparisons must be made on two levels. The differences in methods (differential vs. 
integral) leads to differences in the way geometric terms are handled. These affect 
questions of accuracy and programming efficiency, but are not of a fundamental 
nature. In fact, many algorithms use a combination of both approaches. The 
differences in grids are more fundamental, and affect many problems related to 
computational boundaries. The choice between the two depends on the nature of the 
boundary. Zonal boundaries are more naturally treated with a finite-volume grid. 
In order to achieve strong conservation, a finite-volume grid is also more natural 
for quasi-one-dimensional and axisymmetric flows. On the other hand, greater 
accuracy can be achieved near a topological singularity using a finite-difference 
grid. At a general boundary, such as a solid wall, the choice is not clear cut. Any 
boundary procedure can be adapted to either type of grid. The ultimate choice will 
be determined by programming efficiency and stability considerations. 

There is another class of discretization schemes which utilize both the finite- 
volume and finite-difference grids. Examples of these hybrid schemes are those of 
Ni (Ref. 52) and Roe (Ref. 53) . The conservation law is first applied to the primary 
cells. The changes in conservative variables are then rezoned in a conservative 
manner to yield the changes in the secondary cells. These are then used to update 
the conservative variables in the secondary cells. Since these schemes do not strictly 
fit into either of the classifications according to our definitions, they should properly 
be considered a class onto themselves. 

There is a superficial resemblance between the finite-volume and finite-element 
methods, and much semantic confusion in the literature between the two concepts. 
The author has addressed this question in a previous publication (Ref. 54). Conven- 
tional nodal finite-element methods define the unknowns at cell vertices, and do not 
satisfy the integral conservation laws for those cells. Thus, even though they can be 
formulated in a manner that provides shock capturing, they cannot be related to 
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finite-volume methods. Some recently developed finite-element methods (Refs. 55 
and 56) do satisfy the integral conservation laws. It will be interesting to compare 
their theoretical and practical performance with those of the finite- volume methods. 

APPENDIX: CALCULATION OF ROE AVERAGED EIGENVALUES 

All the calculations in this appendix are carried out for = 0, 7 = 1.4, and 
Pl = PR‘ It follows from Eq. (40) that a = 0.5. 

Case 1. 

In this case the velocity change ur — ul is zissumed to be directed along the 
surface normal n. Given 

cr = 1.2c£, 

u'j^ = 0.995ci (184) 

Ur = 1.205C|;,. 

These conditions define a transonic expansion wave. Using Eqs. (39) and (42), one 
determines for A3 the values 

Ast = — 0.005 c£, 

A 3 R = 0.005CL (185) 

A 3 = -0.005533808ci. 

Note that A3 does not lie between \zl and A3R, 

Case 2. 

By increasing the normal velocity, we can convert Case 1 to 

Cr = \.2cl 

u'i = 1.002CL (186) 

Ur = 1.2102c£,. 

These conditions correspond to a supersonic expansion. The values of A3 are now 

Asl = 0.0002c£, 

A 3 R = 0.0102CL (187) 

A3 = -0.000333808CL. 

Note that the Roe averaged normal velocity is now subsonic, even though both 
and are supersonic. 
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Case 3. 


This case illustrates that anomalous behavior can result just from a velocity 
change tangential to the surface. For concreteness we assume that u^, Ur, and n 
are coplanar. Defining 

ut = lu - Uniil, 

we consider the conditions 

cr = Cl 
u'l, = l.Olcx, 
u'r = 1.02ct 
Wtfi = V-tL + Cl* 

The normal velocity components again define a supersonic expansion. The values 
of A 3 are now 

■^3L = O.Olcx, 

\zR = 0.02ci (tS9) 

A3 = -0.009697516cl. 

The Roe averaged normal velocity for this case is also subsonic. 

These three illustrative calculations indicate that conditions can exist for which 
the Roe averaged eigenvalue lies outside the range determined by the states L and J?. 
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Fig. 3 Boundary half-cell for finite-volume grid. 
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